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ABSTRACT 


The  main  goal  of  this  paper  is  to  discuss  mixed  variational  formulations  for 
time  dependent  problems  such  as  initial  and  boundary  value  problems  for  the  heat 

and  wave  equations  in  a  bounded  domain  O  of  RN(V  >  1).  Then  we  shall 
use  these  formulations  to  derive  mixed  finite  element  approximations  of  the  heat 
and  wave  equations.  Finally^  an  application  to  an  exact  boundary  controllability 
problem  for  the  wave  equation  will  be  presented  together  with  some  numerical 
results.  The  techniques  and  application  briefly  considered  here  will  be  discussed 
with  more  details  in  a  forthcoming  paper. 

INTRODUCTION 

Mixed  variational  principles  and  the  associated  finite  element  approximations 
have  proved  to  be  very  useful  in  order  to  derive  accurate  solution  methods  for 
boundary  value  problems  for  partial  differential  equations.  This  is  particularly  true 
for  elliptic  problems  (see,  e.g.,  [1],  [2]  and  the  references  therein).  A  strong  point 
of  these  techniques  -  compared  to  more  traditional  finite  element  methods  -  is  that 
they  give  fairly  accurate  approximations  of  the  derivatives ;  this  last  property  is  very 
interesting  since  in  many  problems  one  is  more  interested  by  the  derivatives  of  a 
function  than  by  the  function  itself.  Mixed  methods  have  also  been  applied  to  time 
dependent  problems  (see,  e.g.,  [3],  [17])  but  there  are  indeed  very  few  published 
papers  and  applications  where  these  methods  have  been  used  for  time  dependent 
problems  compared  to  the  more  classical  finite  element  methods.  Motivated  by 
optimal  control  applications  (cf.  [4],  [5])  we  shall  briefly  discuss  in  this  short  article 
the  following  topics: 

(i)  Mixed  variational  formulations  for  the  heat  and  wave  equations  (Section  1.). 

(ii)  Mixed  finite  element  approximations  of  the  heat  and  wave  equations  (Section 

(iii)  An  application  to  a  boundary  control  problem  for  the  wave  equation  (Section 
3.). 

1.  MIXED  VARIATIONAL  FORMULATIONS  FOR  THE  HEAT  AND  WAVE  EQUATIONS. 
1.1  Formulation  of  the  basic  time  dependent  problems. 

Let  flbea  bounded  domain  of  R N  (N  >  1)  ;  we  denote  by  T  the  boundary  of 
.  Let  T  be  a  positive  number  (possibly  equal  to  +oo)  ;  we  denote  by  Q  and  E 
the  following  sets  of  Rw+1  : 


Q  =  Q  x  (0,T),E  =  T  x  (0,  T). 

We  suppose  now  that  physical  phenomena  are  taking  place  on  ,  modelled  by 
either  the  following  heat  equation 

(1-1)  ut  —  A u  =  /  in  Q , 

(1-2)  u  =  g  on  E, 

(1-3)  u(x,  0)  =  uQ(x)  on  ft, 
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or  by  the  following  wave  equation 


(1.4)  utt  -  A u  =  f  in  Q, 

(1.5)  u  =  g  on  E, 

(1.6)  u(x,0)  =  u0(x),  ut(x,0)  =  u\(x)  on  Cl. 
In  (1.1)  -  (1.6)  we  have 


x 


r  \N  _ 

{*!' }.•=!»««  —  Qj.  iutt  — 


d2u 

dt*' 


a  =  E 


d2 

dx 2 ' 


It  follows  from,  e.g.  [6],  [7],  that  each  of  the  two  above  problems  has  a  unique 
solution  provided  that  the  data  /  and  a  belong  to  well  chosen  functional  spaces. 
Since  this  paper  is  engineering  oriented  we  shall  not  go  into  the  details  of  those 
(Sobolev  type)  spaces  for  which  the  above  problems  are  well-posed  (there  will  be 
however  some  exceptions). 

1.2  Mixed  variational  formulations  for  problems  (1.1)  —  (1.3)  and  (1.4)  —  (1.6). 

The  key  idea  is  to  take  Vu(V  =  { )  as  master  variable  ;  we  introduce 
therefore  a  new  unknown  p  defined  by 

(1.7)  p  =  Vu(in  Q). 

Assuming  that  u  and  p  are  sufficiently  smooth  we  obtain  -  integrating  by  parts  with 
respect  to  the  space  variables  -  the  following  mixed  variational  formulations: 

Mixed  variational  formulations  of  the  heat  equation  (1.1)  —  (1.3)  : 


(1.8)  /  (ut  —  V  •  p  —  f)vdx  =  0,  VveL2{Cl),a.e.  on  (0,T), 

(1.9)  /  (p  •  q  +  u  V  •  q)dx  =  [  gq  ■  ndT,  WqeH(Q,,  div),  a.e.  on  (0,  T), 

J  n  Jr 

(1.10)  u(x,  0)  =  u0(x)  on  O. 

Mixed  variational  formulations  of  the  wave  equation  (1.4)  —  (1.6)  : 


(1.11)  /  (utt  —  V  •  p  -  f)vdx  =  0,VveL2(£l),a.e.on  (0,T), 

J  Q 

1.12)  j  (p  ■  q  + uW  •  q)dx  =  j  gq  ■  ndT,VqeH(£l,div),a.e.  on  (0,T), 

Jn  Jr 

(1.13)  u(x,  0)  =  u0(x),ut(x,  0)  =  «i(x). 

In  (1.8)  -  (1.13),  we  have  used  the  following  notation:  y  •  z  =  V y,  zelR^;  n 

is  the  unit  vector  of  the  outward  normal  at  T;  dx  =  dx\  •  •  •  dx^  and  finally 

H{Sl,div )  =  {q|qeL2(ft),  V  •  qeL2(H)}. 
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2.  MIXED  FINITE  ELEMENT  APPROXIMATIONS  OF  THE  HEAT  AND  WAVE  EQUATIONS. 
2.1  Generalities. 

With  h  a  space  discretization  step ,  we  approximate  L2(Cl)  and  H(Cl,div )  by  T4 
and  Qh  ,  respectively.  We  suppose  that  14  C  L2(Cl),  Qh  C  H(Cl,  A4)  and  also  that 
Vh  and  Qh  satisfy  compatibility  conditions  implying  convergence  properties  for  the 
corresponding  approximations  (see  e.g.,  [1],  [2]  for  details);  an  important  condition 
to  be  satisfied  is: 

(2-i)  v-QhdVh. 

In  the  particular  case  where  Cl  is  a  2  dimensional  polygonal  whose  boundary  is  the 
union  of  segments  parallel  to  the  coordinate  axis,  we  associate  to  Cl  a  “partition” 

Rh  such  that 


(0 


*-<*>•« -A* 


(ii)  Each  I<  is  a  rectangle  whose  edges  are  parallel  to  the  coordinate  axis, 

(iii)  4f  K  and  A  eRh ,  then  AT  D  K'  =  <f>,  and  either  K  fl  A  —  (f>,  or  K  and  K1  have 
only  a  whole  edge  or  one  vertex  in  common. 

Following  [1],  [2]  and  [8]  -  [10],  a  convergent  choice  for  T4  and  Qh,  constructed 
from  the  above  Rh,  is  given  by: 

(2-2)  Vh  =  {vh\vh\KeQk,VKeRh}, 


(2.3) 


'  Qh  =  {q&|q*  =  {q»ft}f=1,qfcke(Pfc+i  ®  Pfc)  x  (Pk  ®  Pfe+1), 
<  VKeRh]  qih  is  continuous  along  the  edges 
,  of  Rh  parallel  to  0xi+i }; 


m  (2.2),  (2.3),  k  is  a  nonnegative  integer,  Qk  =  Pk  ®  Pk,Pa  is  the  space  of  the 
polynomials  in  one  variable  of  degree  <  s,  and  i  +  1  has  to  be  taken  modulo  2. 
With  such  a  choice  for  Vh  and  Qh  ,  condition  (2.1)  is  clearly  satisfied. 

2.2  Discretization  of  the  heat  equation  (1.1)  —  (1.3). 

Semi  —  Discretization  in  space  : 

Using  the  spaces  Vh  and  Qh  we  shall  “space  discretize”  (1.1)  -  (1.3),  via  (1.8) 
-  (1.10)  as  follows: 

Find  a  pair  {uh(t),ph(t)}eVh  X  Qh,  a.e.  on(0,T),  such  that 


f  du 

(2-4)  ~  ^  *  Pa  ~  fh)vhdx  =  0,VufceI4,  a.e.  on  (0,T), 

(2-5)  qh  +  UhV -qh)dx  =  J  ghqh  •  ndF,\/qheQh,  a.e.  on  (0,T), 

(2.6)  uA(0)=uoft. 


In  (2.4)  -  (2.6),  fh,gh  and  uoh  are  convenient  approximations  of  f,g  and  u0,  respec¬ 
tively  (we  can  take,  for  example,  uQh  as  the  A2 -projection  of  u0  on  Vh)  ■ 
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The  above  approximation  is  not  practical  since  we  still  have  to  solve  an  ordinary 
differential  system,  or  to  be  more  precise  a  system,  coupling  ordinary  differential 
equations  and  (linear)  algebraic  equations. 

Full  Discretization  in  space  —  time  :  Concentrating  (for  simplicity)  on  the  back¬ 
ward  Euler  time  discretization  of  (2.4)  -  (2.6)  we  finally  obtain  the  following  system 
of  difference  -  algebraic  equations  (with  A <(>  0)  a  time  discretization  step): 

For  n  >  0,  find  {u£+1,  p£+1}el4  x  Qh  such  that 

(2.7)  u°h  =  uoh, 

(2'8)  Jn  (“>+A« “  v ' p;+1  “  A”+1)  v>'dx  =  o. 

(2.9)  /  (pj+1  •  q*  +  <+1V  •  q„)  ix  =  f  sjf+'q*  •  ndT,  Vq»eQ*. 

J  n  Jr 

^From  a  practical  point  of  view,  we  can  easily  eliminate  u£+1  from  (2.8),  using 
the  fact  that  V  •  q heVh  ;  we  obtain  then  the  following  linear  variational  equation 
satisfied  by  p£+1  : 

(2  10)  /  -P»+,V-q»+ P»+1  •<»*)<&  =  /r«+1q»-ndr 

l  ~  JnK  +  Ai4”+1)V  •  qhdx^qktQh- p nh+IcQh- 

Solving  (2.10)  can  be  done  by  a  direct  method  -  such  as  Cholesky’s  since  the  bilin¬ 
ear  form  in  (2.10)  is  symmetric  and  positive  definite  -  or  by  a  conjugate  gradient 
algorithm  (see,  for  example,  [11]).  Once  p£+1  is  known,  computing  u£+1  from  (2.8) 
is  straightforward. 

Similarly,  instead  of  backward  Euler,  we  could  have  used  schemes  such  as  for¬ 
ward  Euler,  Crank  -  Nicholson,  multisteps,  Runge  -  Kutta,  .... 

2.3  Discretization  of  the  wave  equation  (1.4)  —  (1.6). 

Starting  from  the  following  variant  of  (2.4)  -  (2.6):  Find  a  pair 

Ph(t)}eVh  x  Qh,  a.e.on(0,T),  such  that 


C2-11)  Ju(~dtF  ~  V  '  Ph  ~  fh)vhdx  -  0,Vufc€Vfc,  a.e.  ora(0,T), 

(2.12)  /  (pfc  •  qA  +  uhV  ■  qh)dx  =  /  ghqh  •  ndT, WqheQh,  a.e.  on(0,T), 

J  a  Jr 

(2.13)  wft(0)  =  uoh,  ^p-(O)  =  uih, 

we  can  fully  discretize  the  wave  problem  (1.4)  -  (1.6)  by  the  following  variant  of  the 
usual  second  order  accurate,  explicit  finite  difference  discretization  scheme  of  the 
wave  equation: 

Assuming  that,  for  n  >  0,u£,p£  and  are  known  compute  first  u£+1  as 

the  solution  of 


(2.14) 


-  V-p£  -  fh)vhdx  -  0yvkeVh;u”+1eVh, 
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and  then  p£+1  as  the  solution  of 


(2.15)  [  p"+1-qtdx=  [  s;+1qi  .  ndT  -  /  uJ+’V  •  qhdx,  Vq„eQh-,  pj+'etj*. 

J  Q  Jt  J£l 

A  most  important  step  is  clearly  the  initialization  of  scheme  (2.14),  (2.15);  assum¬ 
ing  that  f,g,uQ,Ui  are  sufficiently  smooth  we  shall  proceed  as  follows:  compute 

ubuh\Ph  and  ul 

(SAG)  u°h  =  uoh ,  u{  =  uf1 +2Atulh, 


(2.17) 


P  %eQh, 

fa  P°h  '  q hdx  =  fr  g°h qh  ■  ndT  -  fn  uohV  ■  qhdx,VqheQh. 


As  shown  in  [12],  u*(<)  and  ph(t)  will  converge  to  u(t)  and  Vu(t)(u  :  solution  of 
(1-4)  -  (1.6))  as  h  and  At  — »  0  if  a  stability  condition  such  as 

(2.18)  At  <  Ch 


is  satisfied. 

Second  order,  unconditionally  stable  implicit  variants  of  the  above  scheme  can 
be  obtained;  they  will  discussed  in  a  following  paper,  together  with  applications  to 
boundary  control  of  the  wave  equation. 

APPLICATION  TO  AN  EXACT  CONTROLLABILITY  PROBLEM  FOR  THE  WAVE 
EQUATION,  VIA  DIRICHLET  BOUNDARY  CONTROLS. 

3.1  Formulation  of  the  boundary  control  problem. 

We  follow  here  [4],  [5];  we  consider  then  a  phenomenon  taking  place  in  SI  and 
modelled  by  the  wave  equation  (we  keep  the  notation  of  Section  1): 


(3.1) 


utt  —  A u  =  0  in  Q, 


with  the  initial  conditions 


(3-2)  u(x,0)  =  u0(x),ut(x,o)  =  ui(x)  in  SL 

The  problem  here  is  to  find  g  defined  over  £(=  T  x  (0,  T))  such  that  the  following 
final  conditions 


(3-3)  u(x,T )  =  0,  ut(x,T)  =  0  on  S2 

hold  if  one  has 

(3.4)  u  =  g  onE 

as  boundary  condition. 

It  has  been  proved  by  several  authors  (see  [4],  [5],  [13]  for  references)  that  such  a 
g  exists  provided  that  T  is  sufficiently  large  (the  lower  bound  of  the  T’s  for  which 
(3.3)  holds,  Vu0,u:l,  is  -  not  surprisingly  -  of  the  order  of  diameter  (SI)). 
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3-2  Calculation  of  an  exact  Dirichlet  control  via  the  Hilbert  Uniqueness  Method 
of  J .  L.  Lions 

^  [^]>  [5] >  J.L.  Lions  has  introduced  and  analyzed  a  systematic  way  for  con¬ 
structing  Dirichlet  controls  for  which  (3.3)  holds.  The  construction  technique  is 
systematic  and  based  on  the  Hilbert  Uniqueness  Method  (HUM)  to  be  briefly  dis¬ 
cussed  below.  From  now  on,  we  suppose  that 


(3.5) 

u0eL\tt),UleH-\£l)(=  (tf^ft)'), 

where 

aim  -  gv  ££2(si),vi  _  1,.. 

■  N,v  = 

H  *(ft)  is  the  dual  space  of  jEf*(ft), 

and  we 

define  E  and  E'  by 

(3.6) 

E  =  H10(n)xL2(tt),E'  =  H-1(Cl) 

x  L2(ft) 

Next  we  define  an  operator  AeL(E,  E')  as  follows: 

(0 

(it) 

Take  e  =  {e0,ei}e£; 

Integrate  from  0  to  T  : 

(3.7)! 

ftt  -  A<f>  -  0  inQ , 

(3.7)2 

4>  =  0  on 

(3.7), 

^(*,0)  =  eo(x),<f>t(x,0)  =  ex(x) 

on  ft. 

(Hi) 

Integrate  from  T  to  0  : 

(3.8)! 

iftt  —  A  if  =  0  inQ , 

(3.8)2 

(3.8)3 

Tp(x,  T)  =  0 ,ift(x,T)  =  0  on  CL. 
(iv)  take 

(3.9) 

r-C 

O 

1 

cT 

11 

a> 

<J 

where  if(0)  ( resp .  xft( 0)  )  is  the  function  x  — >•  tf(x,0)  ( resp .  x  — >  ipt(x,  0)). 

It  follows  from  J.L.  Lions  [4],  [5]  that  AeL(E,  E'),WT  >  0;  moreover,  if  T  is 
sufficiently  large  (T  >  diameter  (ft)  )  then  A  is  a  strongly  elliptic  operator  from  E 

onto  E  .  In  addition  to  these  properties,  A  is  self-adjoint  and  satisfies  (with  obvious 
notation): 

(3.10)  (A  e.e')=/  f^drA.Ve.e'cE; 

E 

m  (3.10),  (*,  •)  denotes  the  duality  pairing  between  E'  and  E  which  satisfies 

(Ae,e/)=  f(Ae)-e'dx 

Jq 
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if  Ae  is  sufficiently  smooth. 

Application  to  the  exact  boundary  controllability  of  the  wave  equation  : 

(i)  Solve 

(3-11)  Ae  =  {«1,-u0}. 

(ii)  Solve  (3.7),  taking  for  e  ,  in  (3.7)3  ,  the  solution  of  (3.11). 

(iii)  Take  g  =  |J  on  £ . 

If  T  is  sufficiently  large,  it  follows  -  from  the  properties  of  A  -  that  (3.11)  has  a 
unique  solution  in  E  ;  we  have  (cf.  [4],  [5])  and  the  corresponding  solution 

of  (3.8)  satisfies  (3.1)  -  (3.4),  implying  that  g  is  a  Dirichlet  boundary  control  for 
which  the  exact  controllability  property  (3.3)  holds.  Actually,  of  all  the  Dirichlet 
boundary  control  for  which  exact  controllability  holds,  the  one  obtained  by  HUM, 
i.e.  by  solving  (3.11)  is  the  only  one  of  minimal  norm  in  T2(X^))  as  shown  in  [4], 
[5].  From  the  properties  of  A  ,  problem  (3.11)  can  be  solved  by  a  conjugate  gradient 
algorithm  operating  in  space  E  ;  such  an  algorithm  is  described  in  [13],  [14],  together 
with  conforming  finite  finite  element  implementations  of  it. 

3.3  Mixed  formulation  of  the  boundary  control  problem. 

In  fact,  we  shall  describe  a  mixed  formulation  of  problem  (3.11): 

Assuming  that  the  initial  data  u0  and  ui  are  sufficiently  smooth,  so  that  we 
can  use  integral  representations,  the  problem  is  now  to  find  a  triple  {e0,p0,ei} 
satisfying 


(3.12) 


ieo,  PojeWo,  e\eL2(Q,)-,  V{w0, 7r0}eW0,  v\eL2{VL)  we  have 
-  V>(0)ui )dx  =  fn(u !V0  -  u0vi)dx, 


where  in  (3.12): 

(i)  The  space  Wa  is  defined  by 


(3  13)  lW°  =  {N>  7r°}l v0eL2(ty,  n0e(L2(Sl))N,  Jq(tt0  •  q  +  v0  V  •  q)dx  =  0, 
it  can  be  shown  that 

{v0,7 r0}eW0  v0eHl(VL),-K0  =  Vva. 

(ii)  V>(0)  and  ipt(0)  are  obtained  from  e0,p0,ei  as  follows: 


2): 

(3.14), 

[  {4>tt  —  V  •  p)vdx  =  0,  VveL2(0),  a.e. 
Jq 

(3.14)2 

/  (p  •  z  +  •  z)dx  —  0,  Vzei7(0,  div), 

Jq 

(3.14)3 

<f>(x,0)  =  eo(x),^>t(x,0)  =  ei(x)  on  0; 
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then  from  T  to  0  (using  the  fact  that  =  p  •  n  on  ^): 

(3.15) i  j  (tptt  -  V  •  q)vdx  =  0,VueL2(Q),  a.e.  on(0,T), 

JQ 

(3.15) 2  /  (q-z  +  ^V-z )dx  =  I  p  •  n  z  •  ndr,  Vzeif(f2,  div ),  a.e.  on(0,T), 

Jci  Jr 

(3.15) 3  i>(x,  T)  —  0 ,ipt(x,T)  =  0  on  fi.. 

An  easy  calculation  will  show  that  (with  obvious  notation): 

f  JnOM0)eo  “  V»(0)ci  )dx  =  f  p  •  np'  •  ndrdt, 

(3.16)  |  E 

(  V{e0,'K0]ei},{el0,'K,0-,e,1}eW0  x  L2(ft). 

^From  (3.16)  it  appears  that  the  bilinear  form  occuring  in  (3.12)  is  symmetric 
and  positive  semi  definite  ;  actually,  for  T  sufficiently  large  it  is  strongly  elliptic 
(coercive)  over  (W0  x  £2(ft))2.  From  these  properties,  problem  (3.12)  can  be  solved 
by  a  conjugate  gradient  algorithm  operating  in  W0  x  L2(Q,)  ;  such  an  algorithm  is 
described  in  Section  3.4. 

3.4  Conjugate  gradient  solution  of  problem  (S .12). 

3.4.1.  Generalities . 

Problem  (3.12)  is  a  particular  case  of 

(3-17)  Find  ueV  such  thata(u,v)  =  L(v)yveV , 

where  in  (3.17): 

(i)  V  is  an  Hilbert  space,  equipped  with  the  scalar  product  (•,•)  ,  and  the  corre¬ 
sponding  norm  1 1  •  1 1  . 

(ii)  a  :  V  x  V  — >  R  is  bilinear,  continuous  and  V-  elliptic  (i.e.  3a  >  0  such  that 
a(v,v)  >  a||u||2, VueF). 

(iii)  L  :  V  — »•  H  is  linear  and  continuous. 

It  is  well  known  (cf.,  e.g.,  [15,  Appendix  1])  that  under  the  above  hypotheses, 
problem  (3.17)  has  a  unique  solution.  If  in  addition  to  (i)  -  (iii),  the  bilinear  form 

a(v)  is  symmetric  then  problem  (3.17)  is  equivalent  to  the  following  minimization 
one 


(3.18) 


ueV, 

J{u )  <  J(v),VveV, 


with  J(v)  —  2  a(y,u)  —  L(y).  Problem  (3.17),  (3.18)  can  then  be  solved  by  the 
following  conjugate  gradient  algorithm: 

Initialization 


(3.19) 
Solve  then 

(3.20) 


u°eV  is  given. 

g°eV, 

(9°iv)  —  a(u°,v)  —  L(v),WveV. 
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If  9°  —  0,  or  is  “ small ",  take  u  =  u°  ;  if  not,  set 
(3-21)  W°  =  g° ' 

Now  for  n  >  0  ,  suppose  that  un,gn,wn,  are  known  with  wn  ^  0  ;  define  then 
un+1,gn+1,wn+1  as  follows: 

Descent  :  Compute 


(3.22) 

Pn  =  \\gn\\2/< 

and 

(3.23) 

un+1  =  un  - 

Test  of  the  convergence  and  updating  the  descent  direction  :  Solve 


(3.24)  K+W' 

1  (fl,n+1,  w)  =  (gn,v)  -  pna(wn,  v),WveV. 

If  9n+1  =  0  -  or  is  small-take  u  =  un+1  ;  if  not  compute 

(3-25)  7n  =  \\9n+1\\2/\\gn\W 

and  update  wn  by 

(3-26)  wn+1  =  gn+1  +  -jnwn. 

Do  n  =  n  +  1  and  go  to  (3.22). 

The  above  algorithm  converges,  WeV,  and  we  have  (cf.  [16]): 


(3.27) 


|ttn  —  «||  <  Clju0  —  u\ 


m)’’ 


where  C  is  a  constant,  and  where  the  condition  number  ua  is  given  by 


(3.28) 


Va 


sup 

veS 


a{y,v)/ 


inf 

veS 


a(v,v), 


with  S  =  {v\veV,  ||v||  =  1}. 

3-4.2  Application  to  the  solution  of  problem(S .12) 

Since  problem  (3.12)  is  a  particular  problem  (3.17),  with  V  =  W0  x  L2(Q),  it 
can  be  solved  by  the  conjugate  gradient  algorithm  (3.19)  -  (3.26).  An  important 
practical  issue  is  the  proper  choice  of  the  scalar  product  to  be  used  over  W0  x  L2(Q). 
A  fairly  convenient  one  is  provided  by 


(3.29) 


fn(v°vo  +*■<>•<  +  V!  v[)dx, 

V{u0,7r0;  vi},  {v'0,'K'0\v'1}eW0  x  L2(Q). 
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Applying  algorithm  (3.19)  -  (3.26)  to  the  solution  of  problem  (3.12),  with  W0  xL2(fl) 
equipped  with  the  scalar  product  (3.29),  we  obtain  the  following  algorithm: 
Initialization  : 


(3-30)  K,p  DeWo.eleL2^!)  are  given. 

Integrate  then  from  0  to  T  the  wave  equation 

(3-31)i  [  Wt  ~  V  •  p°)vdx  -  0,VueL2(O),  a.e.  on  (0,T), 

(3-31)2  J  (p°  •  z  +  ■  z)dx  =  0,  VzeH(D,div),  a.e.  on  (0,T), 

(3.31)3  f>(0)  =  e°,^(0)  =  e?. 

Then  from  T  to  0  : 

(3-32)i  j  (V’tt  ~  V  •  q°)vdx  =  0,  VveL2(Q),  a.e.  on  (0, T), 

J  Q 

3-32)2  /  (q°  •  z  +  ^ V  •  z)dx  =  [  p°  •  n  z  •  ndr,  Vz eH(Sl,  div), 

Jn  J  p 

a.e.  on  (0,T), 

(3-32)3  V(T)  =  0,W(T)  =  0. 

Compute  then  {g°,ng°}  and  g°  as  follows:  Solve  the  mixed  elliptic  problem: 
Find  {g°0,ng00}eWo  such  that 

(3.33) !  f  (g°0  ~  V  •  ng00)vdx  =  [  (^(0)  -  Ul)vdx,VveL2(Q), 

oft  JQ 

(3.33) 2  j  (7 rg0o  •  q  +  g°S7  ■  q )dx  =  0,  Vqeif(f2,  div), 

J  Q 


and  then 

(3.34) 


g°  =  u0  -  xp°( 0). 


If  {9oi*9o}  {0>0}, g-y  —  0  ,  or  are  small ,  take  p°  •  n|  as  boundary  control;  if 


not,  set 

(3.35) 


E 


=  {g°o,*g0o\g0i}- 


Then  for  n  >  0  ,  assuming  that  {e”,  p">,  ej,  f>n,  ifn,  {g£,  ng?},  g”,  «,  x<},  tof 
are  known,  we  compute  {e£+1,p ”+1},  e?+1,  <^n+1  -0n+1  {an+1  7ran+1  } 
?r+\{<+\™o+1},<+\  as  follows: 

Descent  : 
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Integrate  from  0  to  T 

(3.36) i  [  (<f>u  ~  V  •  P n)vdx  =  0,VveL2(Cl),  a.e.  on  (0 ,T), 

Ja 

(3.36) 2  j  (p"  •  z  +  f  V  •  z)c?x  =  0,  Vzeif(f2,  div),  a.e.  on  (0,  T), 

(3.36) 3  ?”(o)=<,S‘(o)  =  <. 

27ien  /rom  T  to  0  : 

[  (^«  -  V  •  qfWx  =  0,  \/veL2(Q.), 

Jn 

(3.37) i  a.e.  on  (0,  T), 

I  (qn  •  z  +  ^>"V  •  z)dx  =  [  pn  ■  nz  ■  ndT,VzeH (ft,  div), 
Jci  Jr 

(3.37) 2  a.e.  on(0,T), 

(3.37) 3  ?n(T)-0,^n(T)  =  0 

Solve  now  the  mixed  elliptic  problem:  Find  {g'f ,  ng™}eW0  such  that 


(3.38) ! 

(3.38) 2 

and  set 

(3.39) 

Compute  now 

(3.40) 


[  (9o  ~  v  =  [  rfrt(0)vdx,VveL2(£l), 

Jn  Ju 

I  (*9o  ‘  q  +  £0  V  •  q)cte  =  0,  Vqeff(ft,  div), 

Jn 


9i  =  -</>>)• 


=  /n(ig?i2+kg?ia+igria)^ 

fn(w°  &t(0)-v>?W‘(0))dx 

=  /n(lg?l2  +  kgon|2+g"|2)^ 

J  wo  +nJo  -nw?  +Ji  wi)dx 


and  then 

(3.41) 

(3.42) 

(3.43) 

(3.44) 


W+I.P?+1.er+1}  =  {e0”,P;,en-fc{<,  *<,«>?}, 
W‘+i,p”+i}  =  «n,p“}  -MT’.p"}. 

{V”+I,q“+I}  =  {V.”.q"}-p„{r.q"}. 
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Test  of  the  convergence.  New  descent  Direction  : 

V  {9o+1i  7r#o+1>flT+1}  =  {0)0, 0}  -  or  is  small  -  take  pn+1  •  n|  as  boundary 

control;  if  not  compute 


(3.45) 
and  then 


/n(ko+1|2  +  Kff?+1l2  +  \g?+1\2)dx 
fnddo  I2  +  \^9o\2  +  \9i\2)dx 


(3.46)  {w 


n-}-l 

o 


<+'}  =  w 


"+1,>rs.n+1,rf+1 


Do  n  =  n  +  1  and  go  to  (3.36). 

Remark  S .1  :  Problems  (3.33)  and  (3.38)  are  particular  cases  of 


(3.47)i  f  (u  —  V  •  p)vdx  =  f  /udx,  VueL2(ft), 

Jci  J  n 

(3-47)2  I  (p  •  q  +  uV  •  q)dx  =  0,  Vqeff(ft,  div), 

Jn 

which  is  the  mixed  formulation  of  the  following  Dirichlet  problem 
(3-48)  —A u  +  u  —  f  in  0,  u  —  0  on  T. 


Observing  that  V  •  qeL2(0),Vq eH(Sl,div),  we  can  eliminate  u  from  (3.47)i,  (3.47)2 
to  obtain  that  p  satisfies  (if  feL2(tl ))  : 

(3.49)  {PrefJ:n,d™^  ,,  . 

I  Jq(V  '  P^7  •  91  +  P  •  q)<fo  —  —  fn  /V  ■  qdx,Vqeff(Q,  div). 

Solving  (3.49)  (in  fact  its  discrete  variants)  is  fairly  easy  and  can  be  done  by  con¬ 
jugate  gradient  algorithms  (see,  e.g.,  [9]  for  details).  Once  p  is  known,  one  obtains 
easily  u  from  (3.47)!  .  Combining  the  above  algorithm  with  the  mixed  finite  element 
approximations  and  time  discretization  schemes  of  the  wave  equation  discussed  in 

Section  2  is  (almost)  straightforward;  this  issue  will  be  discussed  in  a  forthcoming 
paper. 

3.4.3.  Numerical  experiments. 

The  mixed  finite  element  approximation  and  time  discretization  schemes  of  the 
wave  equation,  described  in  Section  2,  have  been  combined  to  algorithm  (3.30)  - 

(3.46),  to  solve  problem  (3.11)  when  ft  =  (0, 1)  x  (0, 1)  and  T  =  2y/2  .  Using  the 
Fourier  series  techniques  described  in  [13]  we  have  computed  those  initial  data  u0 
and  u i  for  which  the  solution  e(=  {e0,  ei})  of  (3.11)  is  given  by 


(3-50)  e0(x1,x2)  =  sin-Tra;!  sin7rx2,  ei  =  Tv\/2e0. 

We  have  used  the  mixed  finite  element  approximations  of  Section  2,  with  k  —  1 
and  Rh  the  regular  partition  of  ft  associated  to  the  vertices  {ih,jh}  with  0  <  i,j  < 
N,N  being  an  integer  such  that  Nh  =  1  ;  we  have  taken  N  =  16,32,64  .  The 
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time  discretization  of  the  various  wave  equations  involved  in  the  calculations  was 
obtained  using  the  (conditionally  stable)  explicit  scheme  described  in  Section  2. 
Obtaining  the  (approximate)  values  of  the  control  =  p  •  n  on  ^  ,  was  quite 
easy  since  the  values  of  the  fluxes  (i.e.  of  the  normal  components  of  ),  at  the 
element  interfaces  and  at  the  boundary  T  ,  are  the  natural  degrees  of  freedom  for  the 
functions  belonging  to  the  finite  dimensional  space  Qh  approximating  H (Si,  div). 

For  h  =  l/16(resp.l/32, 1/64)  the  finite  dimensional  variant  of  algorithm  (3.30) 

-  (3.46)  converges  in  48  (resp.  72,  119)  iterations  (the  number  of  iterations  varies 

-  approximately  -  like  a /N  ).  These  numbers  are  much  higher  than  those  obtained 
in  [13],  where  the  space  approximation  was  achieved  by  a  conforming  finite  element 
method,  coupled  to  a  biharmonic  TychonofF  regularization  to  eliminate  spurious 
oscillations.  On  the  other  hand,  using,  as  in  the  present  paper,  mixed  finite  element 
approximations,  it  is  not  necessary  to  use  regularization  to  obtain  very  good  nu¬ 
merical  results,  as  shown  in  Figures  3.1  (a),  (b),  (c)  (N=6),  3.2  (a),  (b),  (c)  (N=32), 
3.3(a),  (b),  (c)  (N— 64). 

Figures  (a)  (resp.  (b))  show  the  variation  of  the  exact  (-)  and  computed  (•  •  • )e0 
(resp.  e\  ),  for  0  <  Xi  <  1,X2  =  -5.  Figures  (c)  show  the  variation  on  (0,T)  of  the 
L\V)-  norm  of  the  exact  and  approximate  boundary  controls. 

All  the  above  calculations  have  been  done  on  a  CRAY  X-MP  supercomputer. 

4.  CONCLUSION. 

In  this  paper  we  have  discussed  the  application  of  mixed  finite  element  methods 
to  the  numerical  solution  of  direct  or  inverse  problems  for  time  dependent  equations. 
These  mixed  methods  are  robust  and  accurate.  They  are  however  more  complicated 
to  implement  than  the  traditional  finite  element  methods.  Indeed  many  important 
issues  remain  concerning  the  practical  use  of  the  mixed  methods  considered  here, 
such  as  speeding  up  calculations  by  multigrid  and/or  domain  decomposition  meth¬ 
ods  (cf.  [10]);  we  intend  to  investigate  them  in  the  near  future. 
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Figure  3.1  (c):  Comparison  between  exact  (-)  and  computed  (*)  value  of  ||g(t)||  2 .  (h=l/16). 


Figure  3.3  (b):  Comparison  between  exact  (-)  and 


